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EXECUTIVE  SUMMARY 


A  moving  discontinuous  Galerkin  finite  element  method  with  interface  condition  enforcement  is  for¬ 
mulated  for  flows  with  discontinuous  interfaces.  The  method  preserves  a  high-order  representation  up  to 
the  interface,  without  relying  on  artificial  dissipation  or  shock  capturing  for  stabilization.  The  method  can 
compute  curved  interfaces,  with  a  priori  unknown  topology,  in  multiple  dimensions.  The  method  applies  to 
both  steady  flows,  or  unsteady  flows  via  a  spacetime  formulation. 


E-l 


A  MOVING  DISCONTINUOUS  GALERKIN  FINITE  ELEMENT  METHOD  FOR 

FLOWS  WITH  INTERFACES 


1.  INTRODUCTION 

High-order  methods  have  the  potential  to  resolve  complex  flow  field  structures  with  high  accuracy  and 
at  a  reduced  computational  cost  compared  to  lower-order  methods.  One  such  method  in  particular,  the  dis¬ 
continuous  Galerkin  finite  element  method  (DG)  [1-9],  has  been  used  with  increasing  success  in  the  field 
of  computational  fluid  dynamics  over  the  past  two  decades  to  simulate  a  wide  variety  of  flows  due,  in  no 
small  part,  to  its  ability  to  achieve  high-order  accuracy  on  fully  unstructured  grids.  Another  advantage  of 
the  DG  method  is  a  well-developed  theory  of  adjoint  consistency  [9-12],  which  makes  it  a  powerful  tool 
for  adjoint-based  optimization.  In  addition,  spacetime  discontinuous  Galerkin  methods  [13-15]  provide 
a  unified  discretization  of  unsteady  flow  problems  by  simultaneously  discretizing  space  and  time.  Thus, 
spacetime  DG  offers  the  prospect  of  both  high-order  accuracy  in  space  and  time  and  adjoint  consistency, 
which  provides  an  ideal  framework  for  the  optimization  of  unsteady  flows.  One  question  that  naturally  arises 
is  how  these  remarkable  properties  fare  in  inviscid  flows  that  are  not  smooth  and  contain  discontinuous  in¬ 
terfaces,  including  material  interfaces  and  shocks.  Since  DG  employs  a  discrete  function  space  composed  of 
discontinuous,  piecewise  polynomials,  it  can,  in  principle,  directly  represent  these  discontinuous  interfaces. 
However,  this  requires  that  the  interfaces  are  aligned  with  the  grid.  When  misaligned,  instabilities  can  arise, 
leading  to  the  common  approach  of  artificial  stabilization  or  shock  capturing,  cf.  [16-19].  While  shock 
capturing  has  been  remarkably  successful  for  simulating  a  variety  of  complex  flows,  it  can  diminish  some  of 
the  benefits  of  DG,  including  high-order  accuracy  and  consistency,  an  ongoing  issue  emphasized  in  a  recent 
survey  of  higher-order  methods  by  Wang  et  al.  [20]. 

In  order  to  extend  the  advantages  of  DG  from  smooth  to  piecewise  smooth  flows  with  a  priori  unknown 
interfaces,  we  propose  a  moving  discontinuous  Galerkin  finite  element  method  with  interface  condition 
enforcement  (MDG-ICE).  This  method  fits,  rather  than  captures  interfaces,  so  that  properties  such  as  high- 
order  accuracy  and  adjoint  consistency  might  be  preserved.  There  are  two  important  features  of  MDG-ICE 
that  distinguish  it  from  standard  DG.  First,  the  method  enforces  the  interface  condition  separately  from  the 
conservation  law  in  the  underlying  weak  formulation,  so  that  the  residual  only  vanishes  upon  satisfaction  of 
both.  Second,  the  method  treats  the  discrete  grid  geometry  as  an  additional  solver  variable  so  that  the  solver 
simultaneously  moves  the  discrete  grid  geometry  to  fit  interfaces  while  resolving  the  flow  field.  Thus,  in 
contrast  to  standard  DG,  this  method  has  both  the  means  to  detect,  via  interface  condition  enforcement,  and 
satisfy,  via  grid  movement,  the  conservation  law  and  its  associated  interface  condition. 

1.1  Background 

Interface,  or  more  specifically,  shock  fitting  is  a  technique  with  a  long  history  within  the  field  of  compu¬ 
tational  fluid  dynamics,  as  surveyed  by  Moretti  [21]  and  Salas  [22,  23].  While  a  large  number  of  variations 
exist,  the  basic  idea  of  these  methods  is  to  explicitly  incorporate  the  definition  of  a  shock  surface  into  the 
underlying  grid  structure  and  then  move  the  grid  according  to  the  shock  speed  dictated  by  the  Rankine- 
Hugoniot  conditions.  Recent  work  by  Paciorri,  Bonfiglioli,  et  al.  [24-28]  has  greatly  extended  the  range  of 
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application  of  these  methods  and  clearly  demonstrated  the  advantage  of  preserving  design  accuracy  in  the 
vicinity  of  a  shock,  a  property  sacrificed  by  standard  shock-capturing  methods.  An  ongoing  difficulty  with 
existing  shock-fitting  approaches  is  treating  shocks  or  interfaces  whose  topologies  are  unknown  a  priori  or 
whose  topologies  change  in  time. 

Simulating  material  discontinuities  poses  an  even  more  significant  computational  challenge.  Such  inter¬ 
faces,  which  lack  the  nonlinear  self- steepening  behavior  of  shock  waves,  are  highly  susceptible  to  smearing 
caused  by  artificial  dissipation  in  shock-capturing  methods.  In  an  attempt  to  avoid  this  issue,  a  number  of 
techniques  have  been  developed  that  instead  track  the  motion  of  the  interface  explicitly.  Front  tracking  meth¬ 
ods,  introduced  by  Glimm  et  al.  [29,  30]  and  extended  to  DG  by  Nguyen  et  al.  [31],  analogously  to  shock 
fitting,  are  particularly  well-suited  for  tracking  the  evolution  of  material  or  contact  discontinuities.  Rather 
than  requiring  interfaces  to  be  aligned  with  grid  interfaces,  the  extended  finite  element  method  (XFEM)  [32] 
uses  local  enrichment  functions  to  represent  interfaces  that  require  special  quadrature  rules  to  perform  the 
numerical  integration.  A  broader  family  of  unfitted  finite  element  methods  (FEM)  [15,  33-38]  have  also 
been  developed  for  moving  interface  problems.  An  advantage  of  these  methods  is  that  by  avoiding  the 
alignment  of  grid  interfaces  with  flow  interfaces,  it  becomes  less  onerous  to  simulate  flows  with  unknown 
or  dynamic  interface  topology,  while  providing  high-order  accuracy,  especially  in  comparison  to  interface 
capturing  approaches. 

Lagrangian  or  arbitrary  Lagrangian-Eulerian  (ALE)  methods,  represented  for  higher-order  finite  element 
methods  by  work  such  as  that  of  Dobrev  et  al.  [39-41],  Persson  et  al.  [42],  and  Fidkowski  [43],  are  an¬ 
other  class  of  interface-fitting  techniques.  These  methods  have  numerous  advantages,  particularly  regarding 
highly  accurate  and  robust  treatment  of  dynamic  material  interfaces  or  moving  bodies.  Lagrangian  methods 
excel  at  representing  such  material  interfaces  directly,  as  well  as  moving  bodies,  since  they  are  based  on 
moving  a  spatial  grid  with  the  velocity  of  the  flow.  ALE  methods  further  improve  upon  Lagrangian  meth¬ 
ods,  by  allowing  for  more  flexibility  in  the  grid  velocity.  However,  in  the  presence  of  shocks,  some  form 
of  shock  capturing  appears  to  still  be  required.  In  addition,  for  unsteady  flow  problems,  ALE  methods  can 
deform  spatial  grids  so  severely  that,  between  time  steps,  a  possible  re-interpolation  operation  is  required 
to  preserve  validity  of  the  grid.  In  addition  to  possibly  introducing  error,  such  a  re-interpolation  process 
also  complicates  adjoint  analysis.  In  contrast,  within  a  spacetime  formulation  both  stationary  and  spatially 
deforming  boundaries  are  readily  treated  as  static  spacetime  boundaries,  avoiding  the  difficulties  of  ALE  for 
highly  deforming  geometry  [44],  while  inheriting  the  adjoint  analysis  of  DG  that  has  been  developed  for 
steady  flows. 

Other  related  methods  are  those  that  also  involve  moving  the  grid  geometry  to  improve  accuracy.  Harten 
and  Hyman  [45]  solve  one-dimensional  unsteady  conservation  laws  using  a  moving  grid  based  on  local 
characteristic  velocities.  The  original  moving  finite  element  method  introduces  [46-48]  the  idea  of  treating 
a  spatial  grid  as  a  solver  variable,  which  is  integrated  in  time  for  unsteady  flows.  Roe  and  Nishikawa  also 
proposed  a  fluctuating  splitting  method  which  involves  treating  the  grid  as  a  variable  [49].  Recently,  Sanjaya 
and  Fidkowski  [50]  proposed  an  element  warping  strategy,  which  like  the  present  work  also  moves  element 
nodes  as  part  of  the  solver,  greatly  improving  solution  accuracy,  as  guided  by  least-squares  or  output-based 
error  metrics. 

While  the  aforementioned  interface-fitting  techniques  have  proven  useful  for  a  number  of  flow  scenarios, 
a  number  of  difficulties  still  exist  that  limit  their  general  applicability.  In  particular,  solving  unsteady  flows 
with  complex,  dynamic  interfaces  that  undergo  topological  changes  remains  a  challenge.  The  MDG-ICE 
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method  attempts  to  circumvent  these  difficulties  by  making  interface  fitting  an  intrinsic  part  of  the  underly¬ 
ing  formulation  in  which  the  flow  field  and  interface  locations  are  solved  as  part  of  a  single  unified  method. 
By  utilizing  a  spacetime  formulation,  complex  interface  dynamics  that  can  cause  extreme  grid  distortion 
using  standard  time-marching  interface-fitting  techniques  can  be  handled  by  making  relatively  small  mod¬ 
ifications  to  an  existing  spacetime  mesh.  Furthermore,  incorporating  shock  fitting  directly  into  a  unified 
variational  finite  element  formulation  would  provide  a  natural  setting  to  formulate  adjoint-based  optimiza¬ 
tion  in  the  presence  of  shocks  [51-54],  and  enable  differentiability  with  respect  to  interface  geometry.  Point 
singularities  that  arise  in  inviscid  flows  (e.g.,  rarefactions  and  shock  formations  in  spacetime)  can  also  trig¬ 
ger  instabilities,  due  to  the  inability  of  standard  polynomial  basis  functions  to  represent  the  instantaneous 
transition  that  occurs  at  a  single  point.  While  such  singularities  can  be  stabilized  using  the  same  artificial 
dissipation  introduced  by  shock  capturing,  they  provide  a  source  for  artificial  entropy  generation  which 
contributes  to  the  reduction  in  the  order  of  accuracy  of  the  underlying  numerical  method.  MDG-ICE  also 
provides  a  means  of  naturally  treating  these  types  of  singularities.  Boundary  layers,  which  arise  in  viscous 
flows,  also  require  significant  computational  effort  to  resolve  accurately  and  stably.  Therefore,  a  broader 
goal  beyond  a  method  tailor-designed  for  shocks,  would  be  a  general  method  that  can  exploit  grid  move¬ 
ment  to  help  resolve  a  range  of  challenging  flow  features  including  interfaces,  singularities,  and  boundary 
layers. 

2.  STANDARD  DISCONTINUOUS  GALERKIN  METHOD  FOR  FLOWS  WITH  INTERFACES 

The  discontinuous  Galerkin  method,  which  is  based  on  piecewise  polynomial  approximation  over  po¬ 
tentially  unstructured  grids,  is  an  ideal  candidate  for  computing  piecewise  smooth  flows  with  interfaces. 
As  long  as  any  interfaces  present  in  the  exact  flow  solution  coincide  with  the  interfaces  present  in  the  grid, 
DG  has  the  potential  to  approximate  such  flows  with  high  order  accuracy.  However,  we  seek  a  method  that 
not  only  has  the  means  to  represent  piecewise  smooth  flows,  but  also  provides  the  means  to  detect  inter¬ 
faces.  Therefore,  in  order  to  compute  flows  with  a  priori  unknown  interfaces,  we  seek  a  method  for  which 
the  residual  vanishes  only  upon  satisfaction  of  the  interface  condition,  or  a  weak  or  variational  equivalent 
thereof.  With  such  a  method  at  hand,  the  solver  can  then  introduce  the  discrete  grid  geometry  as  a  variable, 
so  that  as  the  residual  is  driven  to  zero,  the  interfaces  are  fit  with  increasing  accuracy. 

To  demonstrate  the  importance  of  the  residual  vanishing  only  upon  satisfaction  of  the  interface  condi¬ 
tion,  we  consider  a  linear  advection  problem  for  flow  containing  a  jump  discontinuity  in  a  two-dimensional 
domain,  with  an  imposed  uniform  velocity  of  v  =  (vx,  vt)  G  M2.  We  assume  that  vt  =  1  and  interpret  the 
second  dimension  as  representing  time,  while  we  assume  that  the  spatial  velocity  is  positive,  vx  >  0.  There¬ 
fore,  we  can  also  interpret  the  interface  as  an  unsteady,  moving  interface  in  one-dimensional  space.  Since 
the  flow  is  piecewise  smooth,  we  decompose  the  domain  FI  into  two  subdomains  ^1,^2,  separated  by  an 
interface  H  dFl2  with  a  normal  n  =  (nx,nt)  G  M2  that  is  perpendicular  to  the  spacetime  velocity  (vxivt). 
The  two  subdomains  are  defined  as 

Hi  ={(*,*)  <EO.\(x,t)-(nx,nt)  <0} 

&2  =  {(■*, y)  e  (x,t)  ■  (nx,nt)  >0}.  (1) 


The  governing  conservation  law 


V-^-(y)  =  V-(vxy,vty)  =  0,  in£2i,£22 


(2) 
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(a)  The  spacetime  solution  obtained  using  D G(p  =  0).  The  interface  present  in  the  exact  solution  is  illustrated  with  a 
dashed  line. 


(b)  The  piecewise  polynomial  solutions  obtained  for  D G(p  =  0,1,2, 3),  sampled  along  the  line  t  =  1.  The  dashed 
line  located  at  x  =  ^  spans  the  interface  present  in  the  exact  solution. 

Fig.  1:  A  moving  interface  with  spatial  velocity  vx  =  We  observe  that  the  residual  has  vanished  without 
having  fit  the  interface. 
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is  augmented  with  the  initial  condition 


inflow  condition 


y(x,t  =  0) 


2  ^  <  0 

0  jc  >  0  ’ 


y(x  =  -l,t)  =  2, 


(3) 


(4) 


as  well  as  the  interface  condition 


\n  •  (y)]  =  0  on  dFl\  D  <9^2- 


The  exact  solution,  satisfying  Eqs.  (2),  (3),  (4),  (5),  is  the  piecewise  constant  flow 


y(x,t)  = 


2  in  (x,t)  G 

0  in  (xf)  G  ^2 


(5) 


(6) 


We  now  evaluate  the  criterion  that  the  residual  vanishes  only  upon  satisfaction  of  the  interface  condition, 
for  the  case  of  the  standard  discontinuous  Galerkin  method  where  the  grid  is  misaligned  with  the  interface 
present  in  the  exact  solution.  Let  ST  be  a  partition  of  Q,  consisting  of  disjoint  cells  K,  so  that  Q  =  Uke^k. 
Also  let  Y,V  be  piecewise  polynomial  function  spaces  defined  over  cf.  Section  4.1.  The  standard  DG 
weak  formulation  is  given:  find  y  G  Y  such  that 

£  ((h(y+,y~,n),v+)dK-(J?(y),Vv)K)  =0  VveV,  (7) 

where  the  standard  upwind  numerical  flux  for  linear  advection  in 


h(y+,y  ,n) 


(n-v)y~h  if  n  •  v  <  0 
(n-v)y~  if  n  •  v  >  0  ’ 


(8) 


which  is  also  specialized  at  the  domain  boundary  dFl  in  order  to  incorporate  boundary  conditions.  For  a 
detailed  formulation  and  analysis  of  this  method,  for  the  case  of  smooth  flow,  see  Hartmann  and  Leicht  [9, 
11]  and  the  references  contained  therein.  Figure  1  depicts  the  solution  computed  using  DG  for  the  case  of 
vx  =  jq,  along  with  the  intentionally  misaligned  grid  and  interface  present  in  the  exact  solution,  d£l\  n  <9^2- 
Despite  this  misalignment,  the  residual  vanishes  to  machine  precision.  We  observe  that,  in  the  lower-order 
case,  the  solver  produces  a  monotonic  solution,  yet  the  interface  is  diffused  over  many  cells,  which  is 
representative  of  the  stability  of  low-order  methods  even  in  the  presence  of  shocks.  In  the  higher-order  case, 
the  solution  is  less  diffused,  however,  over-shoots  and  under-shoots  are  present,  which  is  indicative  of  the 
instability  of  high-order  methods  in  the  presence  of  discontinuous  interfaces,  necessitating  shock  capturing 
strategies  that  locally  stabilize  the  solution  while,  at  best  locally,  sacrificing  high-order  accuracy  [20]. 
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For  both  high  and  low  order  approximations,  using  the  standard  DG  formulation,  the  residual  vanishes 
without  having  fit  the  interface.  We  therefore  conclude  that  the  standard  DG  formulation  lacks  the  means  to 
unambiguously  detect  a  priori  unknown  interfaces  thus  motivating  the  pursuit  of  a  revised  formulation. 


3.  MOVING  DISCONTINUOUS  GALERKIN  METHOD  WITH  INTERFACE  CONDITION 
ENFORCEMENT 

We  now  present  a  revised  formulation  of  DG,  which  incorporates  the  interface  condition  separately  into 
the  discretization  in  order  to  satisfy  the  criterion  that  the  residual  only  vanishes  once  the  interface  condition 
is  satisfied,  thus  providing  a  means  to  detect  a  priori  unknown  interfaces.  The  method,  which  we  refer  to 
as  the  Moving  Discontinuous  Galerkin  Method  with  Interface  Condition  Enforcement  (MDG-ICE),  is  first 
formulated  in  physical  space  in  Section  3.1.  In  order  to  drive  the  residual  to  zero,  thus  enforcing  both  the 
conservation  law  and  its  associated  interface  condition,  the  method  must  have  a  means  to  fit  interfaces, 
which  is  provided  by  allowing  the  grid  to  move,  or  in  other  words,  treating  the  grid  as  a  variable.  Therefore, 
the  formulation  is  mapped  into  reference  space  in  Section  3.2,  to  clearly  identify  terms  with  a  geometric 
dependence  and  facilitate  differentiation  with  respect  to  geometry. 

3.1  Formulation  in  physical  space  with  fixed  geometry 

3.1.1  Strong  formulation 

Let  Cl  C  W1  be  a  given  domain.  We  assume  that  Cl  is  partitioned  by  ST ,  consisting  of  disjoint  sub- 
domains  or  cells  K*,  so  that  Cl  =  UKe&lc,  with  interfaces  £  composing  a  set  S  so  that  UeG^£  =  U Ke*rdK, 
over  which  a  normal  n  is  defined.  Consider  a  nonlinear  conservation  law,  given  in  strong  form,  defined  for 
piecewise  smooth,  Revalued,  functions  y  as 

V  •  &  (y)  =  0  in  K  Vk*  G 

\n  •  &  (y)]  =  0  on  £  V£  G  S, 

in  terms  of  a  given  flux  function  SE  :  Rm  -A  Rmx”.  In  general  one  might  consider  systems  of  the  form 

v-^(y,Vy)  =  v-(^c(j)-^(yVy))  =^(y),  dD 

where  the  convective  flux  & c  (y)  is  only  a  function  of  the  state  y,  the  diffusive  flux  & d  (y,  Vy)  is  also 
possibly  a  function  of  the  gradient  Vy,  while  5?  (y)  denotes  a  source  term.  In  the  present  setting,  we  restrict 
the  formulation  to  convection  only,  so  that  &  (y)  =  & c  (y)  and  5?  (y)  =0. 

We  assume  that  S  consists  of  two  disjoint  subsets:  the  interior  interfaces 

£q  G  Sq  =  {£q  G  S  |  £q  n  dCl  =  0}  (12) 


(9) 

(10) 


and  exterior  interfaces 


eae^d  =  {ede^\eacdSi} 


(13) 
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so  that  $  =  (Fo  U  .  For  interior  interfaces  £o  £  So  there  exists  /e+,  k  e  ST  such  that  £o  =  <9k*+  H  die  ,  so 
that  (10)  is  defined  as 

0=  =  n+ (y+)+n~  (y~)  :  on  £o  Ve0  €  d?o  (14) 

where  n+,n~  denote  the  outward  facing  normal  of  K*+,  K~  respectively,  so  that  n+  —  —n~.  For  exterior 
interfaces 

0  =  In  ■  &  (y)]  =  n+  ■  &  (y+)  -  rfi  ■  (y+)  ,  on  ed  £d  £  Sd.  (15) 

Here  n+  •  (y+)  is  the  imposed  normal  boundary  flux,  which  may  or  may  not  be  a  function  of  y+  depending 

on  the  type  of  boundary  condition.  Therefore,  we  further  decompose  into  disjoint  subsets  of  inflow  and 
outflow  interfaces  U  <f0 ut,  so  that  at  an  outflow  interface  £out  the  boundary  flux  is  defined  as  the 

interior  flux, 

«+-^out(>’+)  =n+-^(y+),  (16) 

and  therefore  (15)  is  satisfied  trivially,  which  is  equivalent  to  no  boundary  condition  being  enforced.  At  an 
inflow  boundary  £in  £  <fin,  the  normal  boundary  flux  is  a  fixed  constant  independent  of  the  interior  state  y+ 

n+-^-m{y+)=n+-^m.  (17) 


Systems  governed  by  Eqs.  (9),  (10)  include  linear  scalar  advection  (2),  spacetime  Burgers  flow  (69), 
and  inviscid,  compressible  (Euler)  flow,  both  in  steady  and  spacetime  form.  The  Euler  flow  state  variable  is 
given  by 


y=(p,pvh...,pvn,pE)  6Rm, 
where  m  —  n  +  2.  The  i- th  spatial  flux  component  is  given  by 

&i 0)  =  (pv/, pvjV  1  +  p8n pviVn  +  p8jn , pH v, )  €  Mm, 

Here  p  is  density,  (vi , . . . }  v„)  £  M"  is  velocity,  pE  is  stagnation  energy  per  unit  volume,  and 

H  =  (pE  +  p)/p 

is  stagnation  enthalpy.  The  pressure  p  is  defined  as 


(18) 


(19) 


(20) 


(21) 
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where  the  ratio  of  specific  heats  for  air  is  given  as  7  =  1 .4.  The  steady  Euler  flux  is 


while  the  spacetime  Euler  flux  is 

&  (y)  =  (^1  (y) ,  •  •  • ,  (y)  ,y)  g  Mmx(”+1) . 


(22) 


(23) 


3.1.2  Weak  formulation 

Let  the  solution  space  Y  be  a  vector- valued,  broken  Sobolev  space, 


Y  = 


-\yyi  (  r  irn^i 

Hk(^)\  =  {y€  [L2(Q.)]m\VKe^,y\Ke  [Hk  (k)\  }, 


(24) 


where  k  is  a  positive  integer.  The  weak  formulation  is  obtained  by  integrating  the  conservation  law  (9)  and 
its  interface  condition  (10)  against  separate  test  functions:  find  y  €  Y  such  that 


T  (y-^(y)’v)K- L  (["•■^’(y)],w)e  =  0  V(v,w)eVxW,  (25) 

ke£T  EES’ 

or,  based  on  an  integration  by  parts,  equivalently, 

L  ((n-^(y+),v+)^-(^(y),Vv)J-  £([n-^(y)l,w)e  =  °  V(v,w)eVxlT,  (26) 

EES’ 

where  the  test  space  V  =  F,  while  W  is  the  corresponding  single-valued  trace  space. 

In  contrast  to  the  classical  discontinuous  Galerkin  formulation,  we  do  not  substitute  the  flux  function  in 
the  resulting  surface  integral  over  Sk  with  a  numerical  flux,  cf.  [9,  11],  and  instead  retain  the  interior  flux. 
Global  coupling  is  instead  achieved  through  the  weakly  enforced  interface  condition  over  each  £  G  £ .  Since 
the  conservation  law  and  its  interface  condition  are  integrated  against  separate  test  functions,  the  residual 
only  vanishes  upon  weak  satisfaction  of  both.  This  formulation  can  also  be  seen  in  the  analysis  of  adjoint 
flow  problems  involving  shocks  [51-54]. 

3.1.3  Alternative  weak  formulations 

Although  not  considered  in  this  work,  the  weak  formulation  (26)  could  be  modified  to  substitute  a 
numerical  flux,  while  still  enforcing  the  interface  condition  separately:  find  y  G  Y  such  that 

L  {{h{y+W,n)  ,v+)dK-(^(y),Vv)K)  -  £  (ln-^(y)},w)e  =  0  V(v,w)eVx\T,  (27) 

ke£?  EES’ 

which  brings  it  closer  to  the  standard  DG  method.  In  fact,  if  we  choose  IT  to  be  a  zero-dimensional  space, 
then  the  standard  DG  formulation  (7)  is  recovered.  We  observe  also  that  the  choice  of  an  interior  numerical 
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flux  h  (y+  ,  y  ,  n)  =  n-  ffi  (y+ )  recovers  (26),  while  the  exterior  numerical  flux  h  (y+  ,  y  ,n)  =  n-  ffi  (y  )  and 
integrating  by  parts,  would  lead  to:  find  y  G  Y  such  that 

y  ((y-^(y),v)K-  {ln-^(y)j,v+)dK)  -  Y,  (I«-*^(y)l,w)e  =  0  V(v,w)  <EV  xW.  (28) 

ke£T  ee£ 

The  proposed  formulation  is  also  related  to  the  hybridized  DG  formulation  [55,  56],  which  introduces  an 
interface  variable  y  G  W  to  obtain  a  weak  formulation:  find  (y,y)  Gfxf, 

y  ((V-^(y),v)K+  (/i(j+,j,n)  -n-^(j+),v+)aJ-  £  ([A  (y+,J,«)]  ,w)e  =  0  V(v,w)eVxW, 

KE^  EES’ 

(29) 

or 

L  ((^(};+,y,«),v+)aK.-(^'(y),Vv)K)-  £  ([A(y,^,n)l,w)e  =  0  V(v,w)  eV  xW,.  (30) 

ke£?  EES’ 

The  weak  formulation  (26)  is  recovered  by  the  non-standard  choice  of  numerical  flux 

h(y+,y,n)  =n-3?  (y+) .  (31) 

3.1.4  State  equation  and  its  linearization 

We  can  represent  (25)  abstractly  as  a  state  equation  e  (y)  =  0,  by  defining  the  state  operator  e  :  Y  -A 
(V  xW)*foryeFas 

e(y)  =  (v,w)^  Y,  (v-^(y)’v)K-  L  (I”-*^(y)]>w)e>  (32) 

KE &  EES’ 

The  Frechet  derivative  of  the  state  operator  (32)  about  a  reference  state  y  is  then  defined  for  a  given  pertur¬ 
bation  8y  G  Y  as 

e' (y)  8y  =  (v,w)  ^  £  ( V  •  (j)  <5y,  v)  K  -  £  ( [n  •  ^  (y)  5y]  ,w)g.  (33) 

KE^f  EES’ 

3.2  Formulation  in  reference  space  with  variable  geometry 

3. 2. 7  Mapping  from  reference  space 

In  order  to  align  discrete  grid  interfaces  with  flow  interfaces,  the  grid  must  be  treated  as  a  variable. 
Therefore,  we  transform  the  strong  formulation  (9),  (10)  and  weak  formulation  (25)  of  the  flow  equations  in 
reference  coordinates,  in  order  to  facilitate  differentiation  with  respect  to  geometry.  Therefore,  we  assume 
that  there  is  a  continuous,  invertible  mapping 


u  :  FI  — y  £2, 


(34) 
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from  a  reference  domain  Cl  C  W1  to  the  physical  domain  Q  cKn.  We  assume  that  Cl  is  partitioned  by 
so  that  Cl  =  U^e^k.  Also,  we  consider  the  set  of  interfaces  $  consisting  of  disjoint  interfaces  £,  such  that 
=  U ke^dk.  The  mapping  u  is  further  assumed  to  be  (piecewise)  differentiable  with  derivative  or 
Jacobian  matrix  denoted 

Vu\k:  k^3?(Rn,Rn)  Vk*g  (35) 

The  cofactor  matrix  cof  (Vw)  :  k  -A  J2?  (MW,MW) ,  is  defined  for  k  G 

cof  (Vw(x))  =  det  (Vw(x))  (Vw(x))_T  Vx  G  K*.  (36) 

Consider  £  C  $,  for  which  we  assume  that  there  exists  a  parameterization  :  Z)  -a  £,  mapping  from  a 
parameter  space  D  c  A  parameterization  of  £  =  u  (e)  is  then  given  by  the  composition  0£  =  u  o  : 

D  — )►  £.  Given  £  G  ,  the  (non-unit)  normal  s  (Vw)  |g  :  £  — )►  Rn  is  defined  for  x  G  £  as  the  (non-unit)  normal 
of  the  tangent  plane  of  £  corresponding  to  the  parameter  0G1  (x).  For  example,  if  n  =  3,  and  T]  denotes 
the  parametric  coordinate  directions,  then  for  x  G  £ 


5 (Vw)|g  =  6e  A  9e)  o  0g  1 


(37) 


where  d£  A  0£  is  the  vector  (or  cross)  product  of  the  tangent  plane  basis  vectors.  A  general  formula  for 
evaluating  the  vector  product  of  tangent  plane  basis  vectors  is  given  by  the  following:  let  (x\ , . . .  ,xn)  denote 
the  coordinate  directions  in  Rn,  and  the  parameterization  be  given  in  terms  of  components  0£  =  (0e! , . . . ,  0" ) , 
then 


%  0£  A  •  •  •  A  d^_x  0£  =  det 


/  xi 

%  ej 


\  ••• 


(38) 


so  that  in  general 

s  (Vm)  |g  =  (%  0£  A  •  •  •  A  Ge)  o  0A1 , 


(39) 


where  we  assume  also  that  1 1  0$  A  •  •  •  A  1 1  =  1 . 

3.2.2  Geometric  boundary  conditions 

In  order  to  implement  geometric  boundary  conditions,  which  constrain  points  to  the  boundary,  we  intro¬ 
duce  a  nonlinear  projection  operator,  denoted  b:U  ^U.  A  common  and  simple  example  of  such  a  mapping 
would  be  for  the  case  of  constraining  points  to  a  boundary  surface  £  G  <%,  which  lies  on  a  hyperplane 


P  =  {x  G  Rn  \  x-np  =xp  -np}  ,  (40) 

with  a  given  unit  normal  np  and  reference  point  xp.  In  this  case,  we  define  the  projection  to  map  the  point  to 
the  closest  point  on  the  hyperplane 


b(u) |g  —u  —  (np  •  (u—xp))np. 


(41) 
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Another  example  is  the  projection  of  a  point  to  a  fixed  point  jco,  for  example,  the  point  lying  at  the  intersection 
of  n  hyperplanes  in  W \ 

b(u)  |g=x0,  (42) 

so  that  in  particular  the  derivative  of  the  projection  vanishes  at  such  points. 

3.2.3  Strong  and  weak  formulation  in  reference  space 
The  strong  formulation  in  reference  coordinates  is 


(cof  (Vw)  V)  •  (y)  =  0  in  k 

Vice 

(43) 

[5  (Vw)  •  (y)]  =  0  on  £ 

Veei, 

(44) 

b{u)  —  u  —  0  on  e 

V£  E  S 

(45) 

We  assume  that  Y  = 


now  consists  of  functions  defined  in  an  Rm-valued  broken  Sobolev  space 


over  where  k  is  a  positive  integer.  We  also  assume  that  U  =  [H(  (Q)]  ”,  the  Revalued  Sobolev  space  over 
Cl.  We  further  assume  that  the  test  spaces  V  =  Y  and  W ,  the  corresponding  single-valued  trace  space,  now 
consist  of  functions  defined  over  reference  space.  We  then  define  a  provisional  state  operator  e  :Y  xU  -A 

(V  xW)*  for  (y,u)eYxU, 


e(y,u)  =  (v,w)\->  Y,  ((cof(V«)V).^(y),v)j-  £  ([s(Vu) -.^(y)] 


£  ’ 


(46) 


which  has  a  Frechet  derivative  defined  for  perturbation  (8y,  8u)  6f  xU,  and  test  functions  (v,  w)  G  V7  x  W, 
by 


(e' (y,u)(8y,8u),(v,w))  =+  £  (((cof'(Vw)  VSw)  V)  ■  &  (y)  ,v), 

ke^ 

+  L  ((cof(Vw)V)-  (y)  Sy)  ,v)  k 
~  I  ([(^(Vw)V5n)  •  y)]  ,w), 

-  Y :([5(Vm)-  (^'(y)sy)}  (47) 

eei 

The  state  operator  e  :Y  xU  — >  (V  xfff,  which  imposes  geometric  boundary  conditions  (45)  by  composing 
the  provisional  state  operator  (46)  with  the  projection  b  (u),  is  then  defined 

e(y,u)  =  e(y,b(u)),  (48) 

with  Frechet  derivative  defined,  for  state  (y,u)  G  F  x  U  and  perturbation  (Sy,  8u )  G  Y  x  U,  by 


(e'  (y,  u)  (Sy,  8u ) )  =  ey  (y,  b  (u))  8y  +  eu  (y, b  (u))b'  ( u )  Su. 


(49) 
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The  state  equation  in  reference  coordinates  is  e  (y,  u )  =  0.  The  corresponding  weak  formulation  in  refer¬ 
ence  coordinates  is:  find  (y,  u)  G  Y  x  U  such  that 

(e(y,  u) ,  (v,w))  =  0  V(v,w)Gf  xff.  (50) 


4.  SOLVER 

In  this  section,  we  present  an  approach  for  solving  the  MDG-ICE  weak  formulation  in  reference  coordi¬ 
nates  (50),  including  discretization  using  standard  polynomial  finite  elements,  the  nonlinear  solver  strategy, 
and  the  grid  operations  necessary  for  modifying  the  grid  topology  as  the  solver  degenerates  elements. 

4. 1  Discretization 

Discretization  of  the  weak  formulation  (50)  implies  choosing  finite-dimensional  subspaces  Y/7  C  Y,  Uh  C 
[/,  Vh  C  V  and  W/7  C  TV.  We  use  standard  piecewise  polynomials,  cf.  [9],  defined  over  reference  elements. 
Let  S?p  denote  the  space  of  polynomials  spanned  by  the  monomials  xa  with  multi-index  a  G  Nq  ,  satisfying 
Y!i=\  0Ci  <  p.  In  the  case  of  a  simplicial  grid, 

Yh={y€Y\lk€#,y\ft€[&>p]my  (51) 

Let  J3p  denote  the  space  of  polynomials  spanned  by  the  monomials  xa  with  multi-index  a  G  Nq  ,  satisfying 
a{  <  p  for  i  =  1 . . .  n.  In  the  case  of  a  cuboid  grid, 

Yh  =  {yeY\Vke^,y\ke[£p]my  (52) 

The  discrete  subspace  Uh  of  mappings  from  reference  space  to  physical  space  are  also  discretized  into  Re¬ 
valued  piecewise  polynomials,  in  the  case  of  a  simplicial  grid 

Uh  =  G  U  |  Vk  G  w|^.  G  [«^p]  j"  5  (53) 

or  in  the  case  of  a  cuboid  grid 

Uh  =  G  U  |  Vk  G  u\k  G  [<=Sp]  ^  .  (54) 

The  case  that  the  chosen  polynomial  degree  of  Uh  is  equal  to  that  of  Yh  is  referred  to  as  isoparametric.  It  is 

also  possible  to  choose  the  polynomial  degree  of  Uh  to  be  less  (sub-parametric)  or  greater  (super-parametric) 
than  that  of  Yh. 

The  test  space  V/7  is  chosen  to  be  Uh  for  a  Galerkin-type  method.  Wh  is  the  subspace  of  single-valued 
polynomials  defined  over  the  interface.  In  the  case  of  a  simplicial  grid 


Wh  =  (we  w|Vee<#,w|g  e  [^p]m} 


(55) 
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while  in  the  case  of  a  cuboid  grid 


Wh  =  {weW\Vee  £,  w|8  e  [ Qp]m } . 

The  discrete  weak  formulation  is  therefore  just  (50),  only  restricted  to  the  discrete  subspaces, 


(56) 


eh  :  Yh  x  Uh  — ^  Vh  x  W 


(57) 


4.2  Nonlinear  solver 

Since  in  general  dim  (Yh)  +  dim  (Uh)  ^  dim  (Vh)  +  dim  (W/7),  instead  of  solving 

eh  (y,  u)=  0  (58) 


directly,  we  instead  seek  a  stationary  point 


e'h(y,u)  =  0.  (59) 

The  unconstrained  optimization  problem  (59)  can  be  solved  iteratively.  Given  an  initialization  (y,u)0  the 
solution  is  repeatedly  updated 


(y,u)i+1  =  (y,u)i  +  A(y,u)i  i  =  0,1,2,...  (60) 

until  (59)  is  satisfied  to  a  given  tolerance. 

Available  iterative  methods  for  solving  (59)  include  the  first-order  gradient  descent  method  with  incre¬ 
ment  given  by 

S(y,u)  =  -ae'h(y,u)*eh(y,u),  (61) 

the  second-order  Newton  method  with  increment  given  by 

S(y,u)  =  -e'l(y,uyle'h(y,u),  (62) 

or  the  Gauss-Newton  method  with  increment  given  by 

A(y,w)  =  -  {e'h(y,u)*  e'h(y,u))~X  (e'h  (y,u)*eh  (y,u))  .  (63) 

As  the  residual  approaches  zero,  the  Gauss-Newton  method  recovers  the  second  order  convergence  rate  of 
the  Newton  method  (62),  while  avoiding  difficulties  such  as  a  loss  of  positive  definiteness.  Other  choices 
are  possible,  such  as  the  Levenberg-Marquardt  method, 


A  (y,  u)  =  —  (4  (y,  u)*e'h  (y,  u)  +  U  (y,u))  1  (e'h  (y,  u)*  eh  (y,  u)) . 


(64) 


14 


Andrew  Corrigan ,  Andrew  D.  Kercher,  David  A.  Kessler 


where  I  (y,  u )  denotes  the  identity  mapping  and  A  >  0  denotes  a  regularization  factor,  that  ensures  a  full 
rank  linear  system  of  equations.  This  method,  can  be  viewed  as  an  adaptive  combination  of  the  gradient 
descent  (61)  and  Gauss-Newton  (63)  methods,  where  the  former  is  recovered  by  taking  A  large,  while  the 
latter  is  recovered  as  A  — >  0.  Another  possible  choice  would  be  the  Lagrange-Newton  method  for  constrained 
optimization  [57]. 

In  the  present  work,  we  use  the  Levenberg-Marquardt  method  (64),  which  requires  solving  a  symmetric, 
positive  definite  system  of  linear  system  of  equations.  The  regularization  A I  (y,  u)  in  (64)  can  be  generalized 
to  provide  finer-grained  control  of  the  regularization  with  respect  to  different  solution  variables.  Since 
the  system  (59)  is  full  rank  with  respect  to  y,  a  small  or  zero  regularization  can  be  applied  with  respect 
to  y,  partially  recovering  the  Gauss-Newton  method  (63).  On  the  other  hand,  the  system  (59)  is  often  of 
insufficient  rank  with  respect  to  geometry,  so  that  a  greater  regularization  can  be  applied  with  respect  to 
u ,  especially  for  the  initial  nonlinear  solver  iterations.  This  regularization  is  also  useful  for  inhibiting  the 
solver  from  modifying  the  discrete  geometry  u  too  aggressively.  Therefore,  the  regularization  we  use  is  in 
the  present  work  is  of  the  form  (y,  u)  \-A  (Ayy,  A Uu),  with  Ay,  AM  >  0. 

Such  a  positive  definite,  symmetric  system  of  equations  can  be  solved  for  using  either  CG  (conjugate 
gradients)  or  MINRES  [58].  Other  available  least  squares  solvers  include  LSQR  [59]  and  LSMR  [60]. 
The  linear  systems  of  equations  are  in  general  poorly  conditioned,  especially  when  the  regularization  pa¬ 
rameter  is  low,  so  that  preconditioning  is  essential.  In  the  present  work,  we  employ  MINRES  using  right 
preconditioning  with  a  diagonal  preconditioner  [61]. 

4.3  Grid  Operations 

While  in  certain  situations  (cf.  Section  5.1)  the  MDG-ICE  solver  is  able  to  solve  directly  for  the  grid, 
without  modifying  the  grid  topology,  in  general  an  arbitrary  initial  grid  topology  will  not  be  compatible 
with  the  interface  topology  present  in  the  exact  solution  (cf.  Section  5.4).  In  such  situations  MDG-ICE 
drives  cells  in  the  discrete  geometry  u  to  degenerate,  which  can  be  identified  by  their  generating  a  negative 
Jacobian,  detVw  <  0.  In  such  situations,  degenerate  simplicial  cells  can  be  removed  using  standard  local 
grid  operations  including  edge  collapse,  cf.  [62],  while  in  order  to  maintain  a  sufficient  level  of  refinement, 
remaining  cells  can  be  refined  to  maintain  a  target  number  of  cells,  or  to  ensure  that  a  maximum  length 
scale  is  not  exceeded.  In  the  case  of  cuboid  cells,  when  a  cell  degenerates,  we  implement  partial  collapses 
by  aliasing  degrees  of  freedom  on  opposite  faces  of  the  cell,  an  example  of  the  application  of  this  type  of 
degeneration  is  given  in  Section  5.2. 

Upon  convergence  of  the  MDG-ICE  solver  to  a  stationary  point  (59),  if  the  residual  (58)  is  not  yet 
converged  within  a  given  tolerance,  then  the  grid  can  be  refined  and  a  new  stationary  point  sought  on  a  finer 
space  (59). 

4.4  Entropy  Condition  Enforcement 

For  nonlinear  convection  problems,  in  order  to  ensure  uniqueness,  an  entropy  inequality  must  be  en¬ 
forced,  cf.  [63].  In  the  context  of  piecewise  smooth  flow,  this  requires  that 


\n  •  (y)]  >  0  on  £  Ve  G  $ , 


(65) 
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where  her  (y)  denotes  an  entropy  flux.  In  the  case  of  inviscid  compressible  flow,  governed  by  the  Euler 
equations,  the  spatial  entropy  flux  can  be  defined  as 

&s(y)  =  (psvi,...,psv„)  el",  (66) 

while  the  spacetime  entropy  flux  can  be  defined  as 

(y)  =  (psv  1 , . . . ,  psvn,ps )  G  M”+1 ,  (67) 

where  the  entropy  variable,  which  is  negated  in  relation  to  physical  entropy,  can  be  defined  as 

ps  =  -Af-]n(p/p*).  (68) 


In  the  present  work,  a  very  simple  restart  strategy  is  employed  for  the  preliminary  test  cases  considered 
in  Section  4.4.  Between  each  nonlinear  solver  iteration,  along  any  interface  for  which  the  entropy  condition 
is  violated,  the  states  are  averaged.  Since  entropy  condition  violating  solutions  are  dynamically  unstable, 
this  approach  appears  to  be  sufficient  for  the  preliminary  test  cases  considered  in  the  present  work,  as  once  an 
entropy- violating  discontinuity  is  broken,  the  solver  is  not  strongly  driven  to  recreate  it.  For  more  complex 
flows,  a  more  elaborate  entropy  condition  enforcement  may  be  required,  which  will  be  investigated  in  future 
work. 


5.  EXAMPLES 

We  now  use  MDG-ICE  to  compute  a  variety  of  benchmark  flow  configurations,  in  order  to  provide  an 
initial  assessment  of  its  ability  to  stably  and  accurately  compute  piecewise  smooth  flows  with  discontinuous 
interfaces.  The  first  problem,  presented  in  Section  5.1,  a  one-dimensional,  unsteady,  moving  interface,  is 
solved  for  the  case  of  both  a  small  and  large  motion.  This  case  illustrates  the  ability  of  MDG-ICE  to  solve 
for  flows  with  interfaces  with  a  minimal  reliance  on  grid  topology  modification.  While  in  the  first  problem, 
the  interface  is  imposed  as  a  temporal  inflow  condition,  the  second  problem,  presented  in  Section  5.2,  is  that 
of  shock  formation.  This  problem  assesses  the  ability  of  MDG-ICE  to  solve  for  flows  with  grid  topology  that 
changes  in  time.  The  third  problem,  presented  in  Section  5.3,  the  Sod  shock  tube  problem,  assesses  the  abil¬ 
ity  of  MDG-ICE  to  handle  a  number  of  different  types  of  waves  within  a  unified  formulation.  These  include 
a  shock  and  contact  discontinuity,  as  well  as  a  rarefaction,  that  introduces  both  a  singularity  and  derivative 
discontinuities.  The  fourth  problem,  presented  in  Section  5.4,  is  that  of  intersecting  oblique  shocks,  which 
assesses  the  ability  of  MDG-ICE,  in  conjunction  with  local  grid  operations,  to  compute  a  priori  unknown  in¬ 
terfaces  in  both  two-dimensional  and  three-dimensional  space.  The  fifth  problem,  presented  in  Section  5.5, 
of  supersonic  flow  over  a  cylinder,  assesses  the  ability  of  MDG-ICE  to  compute  a  priori  unknown  curved 
interface  geometry,  while  retaining  a  high-order  representation  up  to  the  interface. 

5.1  Moving  Interface 

First,  we  revisit  the  case  involving  a  moving  interface  governed  by  the  spacetime  linear  advection  equa¬ 
tion  with  spatial  velocity  vx  =  ^  considered  in  Section  1.1.  Figure  2  depicts  the  solution  computed  using 
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MDG-ICE.  The  solver  was  initialized  with  a  uniform  grid,  with  the  temporal  initial  condition  (3)  at  t  —  0 
extruded  in  time.  MDG-ICE  computes  the  exact  solution,  and  in  particular,  fits  the  interface  present  in  the 
exact  solution,  which  is  illustrated  with  a  dashed  line.  This  piecewise  constant  flow  was  computed  for  the 
cases  of  p  —  0, 1, 2, 3,  in  order  to  demonstrate  that  MDG-ICE  avoids  the  artificial  dissipation  introduced  by 
low-order  methods  with  unfitted  interfaces,  while  also  avoiding  the  oscillations  produced  by  higher-order 
methods  in  the  presence  of  unfitted  interfaces,  as  shown  in  Figure  1  for  the  case  of  the  standard  DG  method. 
An  analogous  flow  was  computed  using  a  larger  spatial  velocity  vx  =  \  using  MDG-ICE(/?  =  0),  as  shown 
in  Figure  3.  In  this  case,  the  upstream  state  transitioned  across  many  cells,  including  around  one  of  the 
corners  of  the  domain  to  reproduce  the  exact  solution,  without  requiring  any  grid  topology  changes.  This 
case  demonstrates  that  MDG-ICE  is  also  able  to  fit  interfaces  with  large  motion  by  simultaneously  updating 
the  flow  field  and  the  grid,  with  minimal  reliance  on  grid  topology  modification. 

5.2  Spacetime  Burgers  flow 

In  the  previous  example,  an  interface  was  specified  as  an  initial  condition.  An  interesting  question,  is 
whether  or  not  a  method  for  treating  interfaces  can  handle  changes  in  topology  in  time,  which  happens,  for 
example,  in  the  case  of  shock  formation.  An  additional  related  phenomenon  that  occurs  in  inviscid  flows, 
is  that  of  singularities,  which  can  occur  just  prior  to  shock  formation.  To  assess  the  ability  of  MDG-ICE 
to  compute  such  flows,  we  now  consider  a  shock  formation  problem  governed  by  the  spacetime  Burgers 
equation  with  flux 

^(y)=  (^y2,y)  ■  (69) 

The  smooth  initial  condition  at  t  =  0  is  specified  as 

'2  x<\ 

y{x,t  =  0)  =  <  4-8x  \<x<\,  (70) 

,-2 

which  leads  to  a  shock  at  t  =  | .  The  exact  spacetime  solution  is 


2  x  <  x+  ( t ) 

y(.x^)  =  <^Ej  X+ (t)  <x<x~  (t)  , 
,-2  x>x~(t) 

where  the  shock  positions  are  denoted 

x+  (f)  =  minQ+2f,  A  , 


(71) 


(72) 


x  ( t ) 


(73) 
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MDG-ICEQ  =  0) 
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-0.75 


(a)  The  converged  spacetime  solution  obtained  for  =  ^ .  The  solver  was  initialized  with  a  uniform  grid,  with  the 
initial  condition  at  t  =  0  extruded  in  time.  The  solver  was  initialized  with  a  uniform  grid,  with  the  initial  condition  at 
t  =  0  extruded  in  time.  Upon  convergence,  the  interface  present  in  the  exact  solution,  illustrated  with  a  dashed  line, 
has  been  fit. 
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(b)  The  converged  piecewise  polynomial  solutions  obtained  for  MDG-ICE(j?  =  0, 1,2,3),  sampled  along  the  line 
t  =  1.  The  dashed  line  located  at  x  =  ^  spans  the  interface  present  in  the  exact  solution. 


Fig.  2:  A  one-dimensional  moving  interface,  computed  using  MDG-ICE. 
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Fig.  3:  A  one-dimensional  moving  interface,  with  a  velocity  of  vx  =  \  ,  shown  in  a  two-dimensional  space- 
time  domain,  computed  using  MDG-ICE.  The  solver  was  initialized  with  a  uniform  grid,  with  the  initial 
condition  at  t  =  0  extruded  in  time.  Upon  convergence,  the  interface  present  in  the  exact  solution,  illustrated 
with  a  dashed  line,  has  been  fit.  The  interface  exits  the  domain  at  t  —  |. 


This  is  a  variant  of  the  problem  studied  by  Giles  and  Ulbrich  [52-54]  when  considering  adjoint  solutions 
for  shocked  flows. 


The  solution  computed  using  MDG-ICE(j?  =  2)  on  a  grid  discretizing  the  spacetime  domain  Cl  =  (0, 1)  x 
(0,  with  six  quadrilateral  cells  is  shown  in  Figure  4.  The  flow  state  variable  y  was  discretized  using  i?2 
elements,  while  the  shape  variable  u  was  discretized  using  elements.  Figure  4a  shows  the  solution 
initialization,  which  was  the  extrusion  of  the  initial  condition  (70)  in  time.  With  only  six  cells,  MDG-ICE  is 
able  to  reproduce  the  exact  solution  by  partially  degenerating  one  cell,  while  fully  degenerating  another  cell. 
The  degeneration  process  involves  imposing  an  equality  constraint  on  initially  distinct  degrees  of  freedom, 
so  that  upon  subsequent  solver  iterations  they  move  in  a  synchronized  fashion.  The  partially  degenerated 
cell  has  its  top  two  points  collapsed,  enabling  it  to  exactly  reproduce  the  singularity  leading  into  the  point  of 
shock  formation.  As  shown  in  Figure  4c,  the  flow  is  piecewise  linear  in  space,  with  an  increasing  slope  that 
transitions  to  an  instantaneous  jump  at  (jc,t)  =  Q,  |).  Emanating  from  this  point  is  the  fully  degenerated 
cell,  which  has  its  top  two  and  bottom  two  points  collapsed,  coinciding  with  the  shock  interface.  This  case 
indicates  that  MDG-ICE  is  capable  of  handling  interface  topology  that  changes  in  time  as  well  as  flows  with 
singularities. 


5.3  Shock  tube 


In  this  example,  we  consider  a  shock  tube  problem,  governed  by  the  spacetime  Euler  equations  (23), 
where  multiple  interfaces  emanate  from  a  discontinuity  at  t  =  0.  The  initial  conditions  are  those  of  the 
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(a)  Initial  condition  given  by  (70)  extruded  in  time,  with  a 
temporal  grid  spacing  of  At  = 


MDG-ICEfp  =  2) 
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(b)  The  converged  spacetime  solution  computed  by  MDG- 
lCE(p  =  2)  using  a  geometry  representation.  The  point 
of  shock  formation  is  obtained  at  (x,t)  =  ,  |) . 


MDG-ICE(p  =  2) _  _ MDG-ICE(p  =  2) 


(c)  The  converged  piecewise  polynomial  solutions  obtained  for  MDG-ICE  (/?  =  2),  sampled  along  the  lines  t  = 
yg,  |.  The  dashed  line  located  at  (x,t)  =  Q,  g)  spans  the  interface  present  in  the  exact  solution. 


Fig.  4:  Burgers  spacetime  shock  formation 


classic  Sod  shock  tube,  with  a  piecewise  constant  state,  with 


(p,v,p)(x,t  =  0) 


(1,0,1)  x  <  0 
(1,0, J,)  x>0’ 


(74) 


The  solution  computed  using  MDG-ICE(j?  =  3)  on  a  grid  discretizing  the  spacetime  domain  (— x 
(0,  with  ten  quadrilateral  cells  is  shown  in  Figure  5  within  a  single  spacetime  slab.  The  flow  state  variable 
y  was  discretized  using  i?3  elements,  while  the  shape  variable  u  was  discretized  using  elements.  Fig¬ 
ure  5a  shows  the  spacetime  initialization,  which  was  the  extrusion  of  the  initial  condition  (74)  in  time.  Also, 
the  middle  eight  cells  were  degenerated  to  the  point  (x,t)  =  0  at  the  temporal  initialization.  Figures  5c-5f 
show  a  comparison  of  density,  pressure,  velocity,  and  entropy  between  the  MDG-ICE  solution  and  the  exact 
solution,  cf.  [64].  We  observe  in  Figure  5c  that  both  the  shock  and  contact  are  fitted  as  true  discontinuities, 
at  the  correct  location,  indicating  that  the  speed  of  each  is  computed  correctly.  Moreover  at  the  head  and  tail 
of  the  rarefaction,  a  derivative  discontinuity  exists,  which  is  also  fit.  Each  of  these  types  of  discontinuities 
were  computed  without  resorting  to  specialized  logic  specific  to  each  type.  This  highlights  one  of  the  key 
features  of  MDG-ICE,  that  multiple  types  of  discontinuities  can  be  detected  and  fitted  simply  by  solving 
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the  weak  formulation  (50),  which  enforces  both  the  conservation  law  (9)  and  its  interface  condition  (10) 
separately.  Another  key  observation  is  that  even  though  MDG-ICE  does  not  rely  on  artificial  dissipation  for 
stability,  the  computed  solutions  are  oscillation-free  while  in  Figure  5f  we  observe  the  absence  of  spurious 
entropy  generation  at  the  shock. 

5.4  Intersecting  oblique  shocks 


State  Pressure  (p)  Mach  (M)  Density  (p)  Entropy  (s) 

1  I  2  I  1.177652828 

2  1.31540694  1.821253900  1.21557716  1.179751529 

3  1.70550887  1.648699723  1.46260855  1.181536318 

4  2.1 87780929  1 .478144666  1 .74655696  1 . 1 83 1098 10 

Table  1:  Oblique  shocks:  Pressure,  Mach,  Density  and  Entropy  values 
computed  at  each  of  the  distinct  states  in  the  intersecting  oblique  shock 
flow,  shown  in  Figures  6,  7,  and  8  in  2D,  and  Figure  9  in  3D. 


In  this  example,  we  consider  initially  unfitted  intersecting  oblique  shocks,  governed  by  the  steady  Euler 
equations  (22),  in  order  to  assess  the  ability  of  MDG-ICE  to  compute  flows  with  non-trivial  and  a-priori 
unknown  interface  topology.  The  inflow  condition  are  specified  with  a  Mach  number  M\  =  2,  density 
Pi  =  1,  and  pressure  p\  —  1.  The  domain  and  grids  used  for  simulating  this  flow  are  shown  in  Figure  6. 

First  we  computed  this  flow  using  the  standard  DG(p  =1)  method  with  the  HLLC  flux  [64,  65],  aug¬ 
mented  with  a  shock  capturing  strategy  based  on  [16].  The  coarse  grid  with  641  linear  triangular  cells  shown 
in  Figure  6a  was  used  to  compute  the  density  field  shown  in  Figure  7a.  In  order  to  reduce  the  artificial  dis¬ 
sipation  present  in  the  solution,  we  employed  a  residual-based  adaptation  strategy  to  refine  the  grid  locally, 
shown  in  Figure  6b,  resulting  in  reduced  artificial  dissipation  around  the  shocks  as  shown  in  Figure  7b.  We 
also  see  that  the  considerable  amount  of  spurious  entropy  generation  produced  by  shock  capturing  on  the 
coarse  grid,  Figure  8a,  was  reduced  by  adaptive  grid  refinement  as  shown  in  Figure  8b,  but  that  it  is  still 
significant. 

Next  we  computed  this  flow  using  MDG-ICE(p  =  1),  which  does  not  rely  on  any  form  of  shock  capturing 
or  artificial  dissipation.  We  initialized  the  solution  using  the  coarse  grid  DG(p  =1)  shock-captured  solution 
on  641  linear  triangle  cells  (Figure  7a).  As  the  MDG-ICE  solver  progressed,  any  cells  that  degenerated 
were  removed  from  the  grid  via  local  collapse  operations,  as  described  in  Section  4.3.  The  resulting  grid, 
shown  in  Figure  6c,  had  461  linear  triangle  cells  remaining,  with  the  intersecting  shock  interfaces  clearly 
visible.  The  corresponding  density  field,  shown  in  Figure  7c,  shows  that  the  solver  has  fit  the  shocks  as  true 
discontinuities  with  an  instantaneous  jump  across  each  shock  interface.  For  this  piecewise  constant  flow, 
the  MDG-ICE  residual  has  vanished,  indicating  that  the  computed  solution  is  exact.  The  computed  flow 
values,  for  each  flow  field  state  are  listed  in  Table  1.  As  in  the  previous  example,  the  entropy  field,  shown  in 
Figure  8c,  indicates  that  no  spurious  entropy  is  generated. 
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MDG-ICE(p  =  3) 


(a)  The  initial  spacetime  grid  and  density  field  over  Q,  =  (—0.5, 0.5)  x  (0,0.2).  The  initial  solution  consists 
of  10  B3  quadrilateral  cells  using  a  shape  representation,  where  the  middle  eight  cells  are  degenerated 
at  the  point  (x,t)  =  0,  as  part  of  the  temporal  initialization.  The  dashed  lines  indicate  the  location  of  the 
shock,  contact,  and  rarefaction  head  and  tail  that  are  present  in  the  exact  solution. 


MDG-ICE(p  =  3) 


i  ma 


(b)  The  converged  spacetime  grid  and  density  field  over  £1  =  (—0.5, 0.5)  x  (0,0.2).  The  final  solution 
consists  of  10  ^3  quadrilateral  cells  using  a  shape  representation.  The  dashed  lines  indicate  the  location 
of  the  shock,  contact,  and  rarefaction  head  and  tail  that  are  present  in  the  exact  solution. 


(c)  Density  at  t  =  0. 15 


(d)  Pressure  at  t  =  0. 15 

MDG-ICE(p  =  3) 


(e)  Velocity  at  t  =  0.15 


(f)  Entropy  5  =  (y  —  1)  l\n(p/pl)  at  t  =  0.15 


Fig.  5:  Sod  shock  tube 
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Finally  we  computed  three-dimensional  reflecting  planar  shocks  using  MDG-ICE(j?  =  1) ,  with  flow  con¬ 
ditions  that  are  analogous  to  those  of  the  two-dimensional  case.  While,  physically,  the  three-dimensional 
flow  is  equivalent  to  the  two-dimensional  flow,  the  grid  is  fully  unstructured  and  moves  in  all  dimensions. 
The  initial  grid  consisted  of  9,724  linear  tetrahedral  cells.  A  challenge  that  arises  when  computing  three- 
dimensional  fitted  interfaces  is  that  the  grid  topology  has  increased  complexity,  which  can  manifest  in  a 
greater  tendency  for  cells  to  degenerate  and  collapse.  The  final  grid,  which  consisted  of  10,750  linear  tri¬ 
angle  cells,  was  obtained  by  MDG-ICE  in  conjunction  with  local  collapse  operations,  to  remove  degenerate 
cells,  along  with  local  refinement,  to  refine  cells  that  exceeded  a  maximum  length  scale  of  0.1  in  order  to 
ensure  that  sufficient  resolution  was  maintained.  The  computed  flow  field  is  shown  in  Figure  9.  As  in  the 
two-dimensional  case,  the  MDG-ICE  residual  has  vanished,  indicating  that  the  computed  solution  is  exact. 

This  case  shows  that  MDG-ICE  does  not  require  a  priori  knowledge  of  the  shock  location  and  topology, 
or  impose  restrictions  on  the  topology.  As  can  be  seen,  the  method  handles  shock- shock  as  well  as  shock- 
wall  interactions.  This  demonstrates  the  ability  of  MDG-ICE  to  fit  initially  unfitted  shocks,  even  with 
non-trivial  topology.  The  three-dimensional  intersecting  planar  oblique  shocks,  indicates  that  the  ability  of 
MDG-ICE  to  compute  flows  with  non-trivial  interface  topology  generalizes  to  higher  dimensions. 

5.5  Bow  Shock 

So  far  the  problems  considered  have  involved  linear  interface  geometry.  Since  the  MDG-ICE  formula¬ 
tion  makes  no  assumption  about  the  shape  representation,  and  accommodates  standard  finite  element  shape 
representations,  we  consider  a  problem  involving  curved  boundaries  and  interface  geometry.  The  problem 
is  that  of  supersonic  flow  over  a  cylinder,  resulting  in  a  curved  bow  shock  as  described  by  Shu  [66].  Here 
we  also  use  an  initial  grid  consisting  of  15  x  20  points,  or  14  x  19  cells,  which  we  then  decompose  into  532 
unstructured  triangular  cells.  The  inflow  condition  is  Mach  number  Moo  =  3,  pressure  =  1,  and  density 
poo  =  1.4. 

First  we  computed  the  solution  using  the  standard  DG(p  =1)  method  with  the  HLLC  flux  [64,  65], 
augmented  with  a  shock  capturing  strategy  based  on  [16].  This  shock  captured  solution  was  then  used  as  an 
initialization  to  compute  the  flow  field  using  MDG-ICE(j?  =  3)  with  ^-shape  triangle  elements  that  were 
constrained  to  the  outer  ellipse  and  inner  circle  via  geometric  boundary  conditions.  The  grids  are  shown 
in  Figure  10,  along  with  the  computed  pressure  fields  in  Figure  11.  Since  the  exact  solution  has  a  constant 
stagnation  enthalpy  (20)  throughout  the  domain,  Hoo  =  7,  this  quantity  is  shown  as  well  in  Figure  12. 

We  observe  that  the  curved  bow  shock  is  fit  by  the  MDG-ICE  solver  in  conjunction  with  local  grid 
operations  to  collapse  degenerate  cells  while  performing  local  refinement  to  maintain  a  target  number  of 
532  cells.  The  location  of  the  shock  along  the  line  y  =  0  was  computed  as  x  =  —1.698  for  a  stand-off  dis¬ 
tance  of  0.698.  As  a  form  of  shock  fitting  that  does  not  rely  on  any  form  of  shock  capturing  or  artificial 
dissipation,  MDG-ICE  appears  to  be  immune  to  carbuncle  instabilities  [67,  68].  The  computed  pressure 
field  also  appears  to  be  free  of  the  detrimental  effects  of  artificial  dissipation  that  are  apparent  in  the  shock 
captured  solution.  Furthermore,  MDG-ICE  more  closely  preserves  a  constant  value  of  stagnation  enthalpy 
throughout  the  domain,  reducing  the  stagnation  enthalpy  error  \\H  —  Hoo\\L2^  nearly  two  orders  of  magni¬ 
tude  compared  to  the  DG(p  =  1)  solution,  from  1.2633  x  10-1  to  1.9704  x  10-3.  The  pressure  and  Mach 
number  distribution,  sampled  along  the  line  y  =  0,  along  with  the  exact  pressure  at  the  stagnation  point, 
p  «  12.06096,  are  shown  in  Figure  13.  The  computed  pressure  and  Mach  number  at  the  stagnation  point  on 
this  coarse  grid  are  12.0559064  and  2.26753700  x  10-2  respectively. 
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DG(p  =  1)  Coarse  Grid  (641  Cells) 


(a)  An  initial  coarse  grid  of  641  cells. 


(b)  A  finer  grid  obtained  using  isotropic,  residual-based  adaptation,  resulting 
in  25,941  cells. 


(c)  A  fitted  grid  computed  by  MDG-ICE(j?  =1)  without  shock  capturing,  with 
degenerate  cell  collapse,  using  the  coarse  grid  shown  in  (6a)  as  the  initializa¬ 
tion,  resulting  in  461  cells.  The  shock  interfaces  present  in  the  exact  solution 
are  reproduced  by  the  grid. 


Fig.  6:  Intersecting  oblique  shocks:  Grid 
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(a)  The  density  field  computed  using  DG(g  =  1)  with  shock  capturing  on  a 
coarse  grid  with  641  cells.  The  effect  of  artificial  dissipation  introduced  by  the 
shock  capturing  is  noticeable. 


Dftfn  =  n  Rfifinfid  arid  f 25  941  filial 


(b)  The  density  field  computed  using  DG(j?  =  1)  with  shock  capturing  on  the 
finer  grid  with  25,941  cells.  The  effect  of  artificial  dissipation  introduced  by 
the  shock  capturing  is  reduced  from  its  effect  on  the  coarser  grid,  but  is  still 
visible. 


(c)  The  density  field  computed  using  MDG-ICE(/?  =  1)  without  shock  captur¬ 
ing. 


Fig.  7:  Intersecting  oblique  shocks:  Density 


A  Moving  Discontinuous  Galerkin  Finite  Element  Method  for  Flows  with  Interfaces 


25 


(a)  The  entropy  field  computed  using  DG(p  =  1)  with  shock  capturing  on  a 
coarse  grid.  A  significant  amount  of  spurious  entropy  generation  is  visible. 


(b)  The  entropy  field  computed  using  DG(p  =  1)  with  shock  capturing  on  the 
adapted  grid.  Spurious  entropy  generation  has  been  reduced  from  (8a),  how¬ 
ever  is  still  significant. 


(c)  The  entropy  field  computed  using  MDG-ICEQ?  =  1)  without  shock  captur¬ 
ing. 


Fig.  8:  Intersecting  Oblique  Shocks:  Entropy 
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Density 

1.7465- 

1 .4626 


1.2155 

I 


(a)  The  density  field  computed  using  MDG-ICE(p  =  1)  without  shock  capturing:  thresholded  to  show  pi  and  P3. 


(b)  The  density  field  computed  using  MDG-ICE(p  =  1)  without  shock  capturing:  thresholded  to  show  P2  and  P4. 


Fig.  9:  Three-Dimensional  Intersecting  Shocks 
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This  case  indicates  that  MDG-ICE  generalizes  to  the  case  of  curved  interface  geometry.  Moreover,  this 
result  further  illustrates  the  potential  of  MDG-ICE  to  circumvent  the  loss  of  accuracy  that  occurs  using 
shock  capturing. 


6.  CONCLUSIONS  AND  FUTURE  WORK 

The  proposed  method,  MDG-ICE,  accurately  and  stably  computes  flows  with  interfaces,  including 
shocks,  without  relying  on  interface  or  shock  capturing,  within  a  unified  formulation,  by  enforcing  the 
conservation  law  and  its  interface  condition  simultaneously,  while  treating  the  discrete  domain  geometry 
as  a  variable.  Preliminary  test  cases  demonstrate  that  the  method  can  be  used  to  compute  both  steady  and 
unsteady  flows,  a  priori  unknown  interface  topology,  point  singularities,  using  higher-order  elements  in 
arbitrary-dimensional  spaces. 

In  order  to  further  assess  the  robustness  of  the  proposed  method,  future  work  will  apply  MDG-ICE  to 
problems  of  increased  complexity.  For  example,  in  Section  5.1,  a  slowly  moving  discontinuity  was  com¬ 
puted  without  artificial  dissipation  being  introduced.  Analogous  results,  which  are  oscillation-free,  can  be 
obtained  using  MDG-ICE  for  the  case  of  slowly  moving  shocks  in  the  context  of  the  spacetime  Euler  equa¬ 
tions,  which  has  been  reported  as  being  challenging  for  Godunov-type  or  shock  capturing  methods,  cf.  [69- 
71].  Furthermore,  in  such  situations,  unfitted  shocks  can  induce  oscillations  using  standard  methods,  which 
shock  capturing  can  effectively  suppress.  However,  the  challenge  then  becomes  simultaneously  preserving 
other  low-amplitude  flow  structures  such  as  physical  acoustic  waves,  so  that  the  shock  capturing  does  not 
mistakenly  suppress  such  waves  too.  This  issue  is  of  particular  concern  in  aeroacoustics,  and  remains  an 
open  issue  for  higher-order  methods  [20].  This  issue  arises  in  other  applications,  including  in  scramjets, 
in  which  the  flame  can  be  quite  sensitive  to  perturbations  of  the  local  flow  conditions,  or  boundary  layer 
separation  on  hypersonic  vehicles.  In  the  present  work,  we  have  demonstrated  the  ability  of  MDG-ICE  to 
stably  resolve  shocked  flows  without  artificial  dissipation.  In  future  work,  we  will  assess  MDG-ICE’s  ability 
to  resolve  the  interaction  of  smooth  flow  features  with  shocks. 

Future  work  will  further  improve  the  robustness  and  efficiency  of  the  method.  Techniques  such  as  do¬ 
main  decomposition  (for  spacetime  marching)  and  multigrid  are  of  particular  interest  to  improve  the  linear 
solver  efficiency.  The  initialization  and  regularization  strategy  used  by  the  nonlinear  solver  will  be  opti¬ 
mized  to  achieve  robust  convergence  using  a  minimum  number  of  iterations.  Future  work  will  also  consider 
the  extension  of  the  proposed  formulation  to  flows  involving  diffusion  and  reactions,  incompressibility  con¬ 
straints,  as  well  as  other  conservation  laws.  While  in  the  present  work  the  formulation  was  introduced  and 
preliminary  comparisons  demonstrated  the  accuracy  and  stability  of  the  proposed  method,  a  future  theoret¬ 
ical  study  of  convergence  and  stability  will  be  used  to  either  justify  optimal  high-order  convergence  in  the 
presence  of  discontinuous  interface  or  provide  guidance  to  revise  the  method  to  ensure  such  properties. 


Acknowledgements 

This  work  was  sponsored  by  the  Office  of  Naval  Research  through  the  Naval  Research  Laboratory  6.1 
Computational  Physics  Task  Area.  The  authors  acknowledge  Devon  Wood-Thomas  for  help  implementing 
the  grid  refinement  and  edge  collapse  algorithms  and  Aaditya  Singh  for  helpful  discussions  on  linear  solvers 
and  preconditioners.  In  addition,  the  authors  acknowledge  Rainald  Lohner,  Johann  Dahm,  Brian  Taylor, 
Ryan  Johnson,  Doug  Schwer,  Ravi  Ramamurti,  David  Mott,  and  Gopal  Patnaik  for  helpful  discussions  and 
providing  feedback  on  this  work. 


28 


Andrew  Corrigan,  Andrew  D.  Kercher,  David  A.  Kessler 


DG(v  =  n 


X 

(a)  532  triangle  cells 


MDG-ICE(p  =  3) 
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x 

(b)  532  curved  ^2  triangle  cells 


Fig.  10:  The  initial  linear  and  quadratic  grids,  as  well  as  the  final  shock  fit  grid,  for  the  bow  shock  simulation. 
The  DG(/?  =  1)  solver  used  shock  capturing  on  a  grid  with  532  triangle  cells.  The  MDG-ICE(/?  =  3) 
solver  was  initialized  by  projecting  the  DG(j?  =1)  solution,  and,  in  conjunction  with  degenerate  cell  collapse 
and  local  refinement,  successfully  fit  the  shock  using  curved  ^2  triangle  cells,  without  shock  capturing. 
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MDG-ICE(p   3) 
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(a)  532  ?A\  triangle  cells  (b)  532  -field/ -shape  triangle  cells 


Fig.  11:  The  computed  pressure  field  corresponding  to  the  grids  shown  in  Figure  10.  The  MDG-ICE(p  =  3) 
solver  was  initialized  by  projecting  the  DG(p  =  1)  solution.  The  flow  field  variables  were  resolved  using  ^3 
elements,  without  shock  capturing,  while  the  shock  interface  was  fit  using  curved  a  .  z’i  shape  representation. 
The  location  of  the  shock  along  the  line  y  =  0  was  computed  as  x  =  —  1 .698  for  a  stand-off  distance  of  0.698. 
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(a)  532  triangle  cells  (b)  532  ^3 -field/ ^2 -shape  triangle  cells 


Fig.  12:  The  computed  stagnation  enthalpy  (20)  field  corresponding  to  the  grids  shown  in  Figure  10.  The 
exact  solution  has  a  constant  stagnation  enthalpy  Hoo  =  7.  The  stagnation  enthalpy  field  computed  using 
MDG-ICE(/7  =  3)  without  shock  capturing,  significantly  improves  the  accuracy  with  respect  to  the  exact 
solution,  decreasing  the  error  \\H  —  fto||L2(^)  decreased  from  1.2633  x  10-1  to  1.9704  x  10-3. 
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MDG-ICE(p  _=  3) 


(a)  The  pressure  sampled  along  y  =  0.  The  location  of 
the  shock  was  computed  as  x=  —1.698  for  a  stand-off  dis¬ 
tance  of  0.698.  The  exact  pressure  at  the  stagnation  point, 
p  «  12.06096,  is  marked  with  a  cross.  The  computed  pres¬ 
sure  at  the  stagnation  point  on  this  coarse  grid  is  12.0559064. 


(b)  The  Mach  number  sampled  along  y  =  0.  The  exact  Mach 
number  at  the  stagnation  point,  M  =  0,  is  marked  with  a  cross. 
The  computed  Mach  number  at  the  stagnation  point  on  this 
coarse  grid  is  2.26753700  x  10-2. 


Fig.  13:  Bow  shock:  Computed  with  MDG-ICE(/?  =  3)  using  a  ^2  shape  representation. 


A  version  of  this  work  has  been  submitted  for  consideration  for  publication  in  the  Journal  of  Computa¬ 
tional  Physics. 
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